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ABSTRACT 

We consider the interface between an accretion disk and a magnet osphere 
surrounding the accreting mass. We argue that such an interface can occur not 
only with a magnetized neutron star but also sometimes with an unmagnetized 
neutron star or a black hole. The gas at the magnetospheric interface is generally 
Rayleigh- Taylor unstable and may also be Kelvin-Helmholtz unstable. Because 
of these instabilities, modes with low azimuthal wavenumbers m are expected 
to grow to large amplitude. It is proposed that the resulting nonaxisymmetric 
structures contribute to the high frequency quasi-periodic oscillations that have 
been seen in neutron-star and black-hole X-ray binaries. The mode oscillation 
frequencies are calculated to be approximately equal to mQ m , where fl m is the 
angular velocity of the accreting gas at the magnetospheric radius. Thus, mode 
frequencies should often be in the approximate ratio 1:2:3, etc. If the pressure of 
the gas in the disk is not large, then the m = 1 mode will be stable. In this case, 
the mode frequencies should be in the approximate ratio 2:3, etc. There is some 
observational evidence for such simple frequency ratios. 

Subject headings: accretion, accretion disks — black hole physics — magnetic 
fields — MHD — instabilities — stars: oscillations — X-rays: binaries — X- 
rays: stars 

1. Introduction 

The study of the timing properties of X-ray binaries and, in particular, quasi-periodic 
oscillations (QPOs) in these sources, has for long been a major area of research (see van der 
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Klis 2000 and Remillard et al. 2002b for recent reviews, and Lewin, van Paradijs, & van der 
Klis 1988 for a discussion of earlier observations). The field has become especially active in 
recent years after the launch of the Rossi X-Ray Timing Explorer (RXTE; Bradt, Rothschild, 
& Swank 1993) with its superb timing capability. Among the divergent phenomenology of 
QPOs observed in X-ray binaries, of particular interest are the "kilohertz QPOs," which 
have frequencies ~ 10 2 — 10 3 Hz. 

About twenty accreting neutron-star sources and five accreting black-hole sources are 
known to exhibit kHz QPOs (van der Klis 2000; Remillard et al. 2002b). These QPOs often 
appear in pairs (so called "twin QPOs"), with the frequency difference being anticorrelated 
to the QPO frequencies (van der Klis 2000; Di Salvo, Mendez, & van der Klis 2003, see, 
however, Migliari, van der Klis, & Fender 2003). In many neutron star systems the frequency 
difference appears to be related to the spin frequency of the neutron star (van der Klis 2000; 
Wijnands et al. 2003). At the same time, in several sources, the twin kHz QPO frequencies 
seem to occur in the ratio of simple integers, with a frequency ratio of 2:3 being common in 
both black hole and neutron star systems (Abramowicz & Kluzniak 2001; Abramowicz et al. 
2003; Remillard et al. 2002a). 

Kilohertz QPOs occur at a similar frequency in neutron star sources that differ in X- 
ray luminosity by more than two orders of magnitude. Moreover, the QPO frequency in 
a given source seems to track the fluctuation in the X-ray intensity from average for that 
source rather than the absolute luminosity of the source (van der Klis 1997, 1998, 2000; 
Zhang, Strohmayer, & Swank 1997). Furthermore, there is a surprising continuity of QPO 
properties between black hole and neutron star binaries (Psaltis, Belloni & van der Klis 1999; 
Belloni, Psaltis & van der Klis 2002). The continuity even appears to extend to white dwarfs 
(Warner, Woudt & Pretorius 2003). 

Because of their high frequencies, kHz QPOs must be produced by processes close to 
the accreting mass. However, since these QPOs have been observed in both neutron-star and 
black-hole X-ray binaries, the oscillations are unlikely to be associated with the surface of the 
accreting object. Instead, it seems likely that the kilohertz QPOs originate in the accretion 
flow surrounding the central mass. Motivated by this argument, a variety of accretion-based 
models have been proposed to explain the oscillations (Lewin, van Paradijs, & van der Klis 
1988; van der Klis 2000, and references therein). 

The classical beat frequency model for QPOs in neutron-star sources proposes that one 
of the kHz QPOs is associated with the Keplerian frequency of the orbiting gas at the inner 
edge of the accretion disk and that the second kHz QPO represents a beat phenomenon 
between the orbiting gas and the spin of the central star (Alpar & Shaham 1985; Lamb et al. 
1985; Strohmayer et al. 1996; Miller, Lamb, & Psaltis 1998; Lamb & Miller 2001). This model 
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predicts that the frequency difference between the two kHz QPOs should be equal to the 
spin frequency of the neutron star and should be constant. However, observations show that 
the separation is not a constant: it generally decreases when the QPO frequencies increase. 
The frequency difference is also observed not to be precisely equal to the frequency of burst 
oscillations (van der Klis 2000), which are believed to match the spin frequency of the star. 
Most interestingly, kHz QPOs have recently been detected in an accreting millisecond pulsar, 
SAX J1808.4-3658, whose spin frequency is known (Wijnands et al. 2003). The ~ 195-Hz 
frequency difference between the two kHz QPOs in this source is far below the 401-Hz spin 
frequency of the neutron star, but is consistent with half the spin frequency. Wijnands et 
al. (2003) argue that their observations falsify the beat-frequency model and pose a severe 
challenge for all current models of the kHz QPOs (but see Lamb & Miller 2003). 

Stella & Vietri (1998,1999) proposed that the QPOs may be interpreted in terms of 
fundamental frequencies of test particles in motion in the relativistic potential of a neutron 
star or a black hole. They identified individual QPO frequencies with the orbital, periastron 
precession and nodal precession frequencies of test particles and showed that some predicted 
scalings between the frequencies agree with observations (Stella & Vietri 1999; Stella, Vietri 
& Morsink 1999). In contrast to some versions of the beat frequency model, this model does 
not require a magnetic field anchored in the central star, and thus provides a natural expla- 
nation for why twin QPOs are seen in both neutron-star and black-hole systems. However, 
the fact that the frequency difference between twin kHz QPOs in neutron star systems is 
often of order the neutron star frequency has no simple explanation. Also, since the model 
explicitly invokes general relativistic effects, it cannot explain the continuity in QPO prop- 
erties between neutron stars and black holes on the one hand and white dwarfs on the other 
(Warner et al. 2003). 

A number of scientists have investigated modes of oscillation of an accretion disk as a 
model of high frequency QPOs (for reviews see Kato, Fukue, & Mineshige 1998; Wagoner 
1999; Kato 2001a). In an important paper, Okazaki, Kato, & Fukue (1987) showed that 
axisymmetric g-modes are trapped in the inner regions of a relativistic disk, where the 
epicyclic frequency k reaches a maximum. This idea was exploited by Nowak & Wagoner 
(1991, 1992) and a number of other workers (Perez et al. 1997; Nowak et al. 1997; Silbergleit 
et al. 2001; Wagoner, Silbergleit, & Ortega- Rodriguez 2001; Abramowicz & Kluzniak 2001) 
who worked out the physical properties of the trapped modes. Much of the work has focused 
on neutral (i.e., non-growing) modes. The idea in such models is that, although the modes 
are stable, they might nevertheless be excited by some driving mechanism such as disk 
turbulence to produce the observed QPOs. 

Recently, Kato (2001b, 2002) claimed to find that nonaxisymmetric g-modes are trapped 
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between two forbidden zones that lie on either side of the corotation radius and that the 
modes are highly unstable. Such dynamically unstable modes are very interesting since 
they would be self-excited and would spontaneously grow to non-linear amplitudes without 
the need for an external driving mechanism. Motivated by Kato's work, Li, Goodman, & 
Narayan (2003) studied nonaxisymmetric g-mode and p-mode instabilities in an unmagne- 
tized isothermal accretion disk. They found that g- and p-modes with a nonzero number of 
vertical nodes are strongly absorbed at corotation and thus not amplified. Waves without 
vertical nodes are amplified (as known earlier, see Papaloizou & Pringle 1985; Goldreich, 
Goodman, & Narayan 1986; Narayan, Goldreich, & Goodman 1987), but the amplification 
is very weak. Therefore, any energy loss, either during propagation of the wave or during 
reflection at the boundaries, would kill the instability. Thus, Li et al. (2003) concluded, in 
agreement with Kato (2003), that nonaxisymmetric disk modes are not sufficiently unstable 
to be of interest for the QPO problem. Ortega-Rodriguez & Wagoner (2000) found that 
viscosity causes the fundamental g- and p-modes in a disk to grow, though the relevance of 
their result for QPOs is presently unclear. 

Since the amplitude of intensity fluctuations observed in QPOs is often quite large, it 
is the opinion of the present authors that the oscillations are likely to be the result of a 
strongly-growing instability of some sort in the accretion flow. Having argued above that a 
hydrodynamic instability is ruled out (Kato 2003; Li et al. 2003), it is natural to consider the 
effect of magnetic fields via a magnetohydrodynamic (MHD) instability. Since, in general, 
only a small number of frequencies dominate the intensity fluctuations, the instability cannot 
be present at all disk radii (as in the case of the magnetorotational instability, Balbus & 
Hawley 1998), but must be associated with some special radius in the disk. The most 
natural choice for the special radius is the inner edge of the disk. Motivated by these 
arguments, we consider in this paper an accretion disk that is terminated at its inner edge 
by a strong vertical magnetic field. The resulting magnetospheric interface suffers from both 
the Rayleigh- Taylor and Kelvin-Helmholtz instabilities. We calculate the frequencies and 
growth rates of the unstable modes and consider the role that these modes play in the kHz 
QPOs seen in X-ray binaries. 

The Rayleigh- Taylor instability and the related interchange instability have been con- 
sidered by a number of authors, both for modeling QPOs (Titarchuk 2002, 2003) and in 
connection with other applications (Arons & Lea 1976a,b; Eisner & Lamb 1976, 1977; Michel 
1977a,b; Spruit & Taam 1990; Kaisig, Tajima & Lovelace 1992; Lubow & Spruit 1995; Spruit, 
Stehle & Papaloizou 1995; Chandran 2001). We compare our work to these earlier studies. 

The plan of the paper is as follows. In §2 we describe the model and its initial equi- 
librium state. In §3 we consider linear perturbations, derive a wave equation satisfied by 
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the perturbations, and obtain the corresponding jump conditions across the magnetospheric 
radius where there is a discontinuity in disk properties. In §4 we consider the case when 
the mass density of the disk has a jump but the angular velocity is continuous across the 
boundary; this situation occurs when the central object is a black hole. We show that, under 
suitable conditions, there is a Rayleigh- Taylor instability at the boundary, and we study in 
detail the growth of the instability using both analytical and numerical methods. In §5 we 
study the case when both the mass density and the angular velocity of the disk are discon- 
tinuous at the boundary, which corresponds to the situation when the central object is a 
neutron star. We show that both the Rayleigh- Taylor and the Kelvin-Helmholtz instabili- 
ties are present under suitable conditions. In §6 we summarize and discuss the results, and 
briefly compare the predictions of the model to observations of QPOs. In Appendix A we 
derive a necessary condition for disk instability. In Appendix B we briefly review classical 
results on the Rayleigh- Taylor and Kelvin-Helmholtz instabilities in Cartesian flows. 

2. Basic Model and Initial Equilibrium 

2.1. Outline of the Model 

We consider a differentially-rotating axisymmetric accretion flow that is terminated on 
the inside by a strong magnetic field. The interface between the disk and the magnetosphere 
occurs at the magnetospheric radius r m . The density experiences a jump at r m , being larger 
on the outside and smaller on the inside. There may also be a jump in the angular velocity. 
Because of the density and velocity jumps, the system is potentially Rayleigh- Taylor and 
Kelvin-Helmholtz unstable (Chandrasekhar 1961; Drazin & Reid 1981). 

The geometry we envisage is natural for accretion onto a magnetized neutron star. This 
geometry has been studied for many years (e.g., Ghosh & Lamb 1978, 1979) and there are 
several discussions of possible instabilities at the disk-magnetosphere interface (Arons & Lea 
1976a,b; Eisner k Lamb 1976, 1977; Michel 1977a,b; Ikhsanov & Pustil'nik 1996). It has 
been suggested that the interchange instability (i.e., the Rayleigh- Taylor instability in the 
presence of magnetic fields) is the fastest mode of mass transport into the magnetosphere 
of an accreting magnetized neutron star (Eisner & Lamb 1984). When the accretion flow 
is spherical and the neutron star rotates slowly, the instability criterion is determined by a 
competition between the gravitational force of the neutron star and the poloidal curvature 
force of the magnetic field: gravity tends to drive the Rayleigh- Taylor instability, while the 
curvature force tends to stabilize the configuration (Arons & Lea 1976a,b; Eisner & Lamb 
1976, 1977; Michel 1977a,b). When the neutron star rotates fast enough, the magneto- 
spheric boundary becomes stable with respect to the interchange instability because of the 
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strong toroidal magnetic field that is produced in the boundary layer by the rotation of the 
magnetosphere relative to the spherical accretion flow (Ikhsanov & Pustil'nik 1996). 

Similar processes for a thin Keplerian disk accreting onto a magnetized neutron star 
have been analysed by Spruit & Taam (1990). These authors confirmed that the Rayleigh- 
Taylor instability can be very effective in transporting mass across field lines, allowing the 
gas to drift along the midplane in a disk-like fashion even in regions of the magnetosphere 
that corotate with the central star. The interchange instability in a general thin accretion 
disk with differential rotation has been studied by Lubow & Spruit (1995) and Spruit et al. 
(1995). They found that disk shear tends to stabilize the configuration. 

The model we consider differs from the studies cited above in one important respect: we 
explicitly avoid field curvature and focus on the pure Rayleigh- Taylor and Kelvin-Helmholtz 
problem for a disk-magnetosphere interface. The model is thus highly simplified and is per- 
haps most relevant for gas in the mid-plane of the accretion disk. We also ignore the toroidal 
component of the magnetic field, even though it may often be important (Ghosh & Lamb 
1978; Ikhsanov & Pustil'nik 1996). A toroidal magnetic field affects the stability of the disk 
in two ways: (1) The toroidal field causes an inward force because of its curvature, enhancing 
the effective gravity in the disk and thereby boosting the Rayleigh- Taylor instability. (2) The 
toroidal magnetic field also produces a tension at the interface which tends to stabilize the 
disk (see, e.g, Chandrasekhar 1961). The overall effect will be determined by the competition 
between these two tendencies. We expect that for modes with a low azimuthal wave number 
the two effects will approximately cancel each other so that, at least for modest levels of 
toroidal field, the results may not be affected significantly. However, for modes with a large 
azimuthal wave number, the second effect will be dominant and so a toroidal field will very 
likely stabilize the disk. 

Although our model is clearly relevant for neutron star accretion, it is in fact relevant also 
for accretion onto a black hole or an unmagnetized neutron star. Bisnovatyi-Kogan & Ruz- 
maikin (1974, 1976) proposed that an accretion flow around a black hole may, under certain 
circumstances, advect a considerable amount of magnetic flux to the center. Once enough 
magnetic flux has collected, the accumulated field will disrupt the flow at a magnetospheric 
radius, exactly as in our model. This idea has been confirmed in recent three-dimensional 
MHD simulations of radiatively inefficient accretion flows by Igumenshchev, Narayan, & 
Abramowicz (2003), and some additional consequences have been explored by Narayan, Igu- 
menshchev & Abramowicz (2003). Something like a magnetosphere may also result from a 
global poloidal magnetic field being generated within the disk (Livio, Pringle, & King 2003). 

The angular velocity profile has some important qualitative differences between the 
black hole and neutron star problems. In the black hole case, we expect the angular velocity 
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to be continuous across r m (§4, Fig. 2), whereas in the neutron star case the angular velocity 
will in general be discontinuous at r m (§5, Fig. 7). Thus, the disk-magnetosphere interface in 
the case of a black hole experiences only the Rayleigh- Taylor instability, whereas in the case 
of a neutron star it can have both the Rayleigh- Taylor and Kelvin-Helmholtz instabilities. 

In order to be able to handle both the black hole and neutron star problems, we set up the 
theoretical framework sufficiently generally as far as the angular velocity profile is concerned. 
However, to make the analysis tractable, we make some other drastic simplifications. First, 
we assume that the gas is incompressible. With this approximation, compressive waves are 
filtered out. The interchange and Kelvin-Helmholtz instabilities, which are our main interest 
and which are both incompressible in nature, survive unaffected (Spruit et al. 1995). Second, 
we assume that the magnetic field is purely along the z-axis (i.e., parallel to the rotation 
axis of the gas flow) and that all fluid quantities are independent of z. We thus solve a 
cylindrically symmetric problem in which d/dz — and all quantities are function only of 
the cylindrical radius r and the azimuthal angle 0. The geometry of the model is as shown 
in Figure 1. 

We describe below in §2.2 the equilibrium flow, and analyse in §3 the effect of linear 
perturbations. 



2.2. Basic Equations and the Initial Equilibrium State 

The dynamics of a perfectly conducting and magnetized Newtonian accretion disk is 
governed by the basic equations of MHD (Balbus & Hawley 1998, we neglect viscosity and 
resistivity): 

^ + V.(pv)=0, (1) 

"(! +v ' v ) v =-" v *- v ( p+ £) + il B ' vB - < 2 > 

— = V x (v x B) , (3) 
V • B = , (4) 

where p is the mass density, v is the velocity, $ is the gravitational potential of the central 
compact object, p is the gas pressure, and B is the magnetic field. Equations (l)-(3) are 
respectively the continuity, momentum and induction equations, and equation (4) is one of 
Maxwell equations. 
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We assume that the disk fluid is incompressible, i.e. 

V • v = . (5) 

Substituting this in the continuity equation (1) gives 

dp d d „ 

where d/dt is the Lagrangian time derivative as measured in the comoving frame of the fluid. 
Equation (6) states that the mass density measured by a comoving observer does not change 
with time, which is obvious for an incompressible fluid. 

As explained in §2.1, the equilibrium state is assumed to be described by 

P = Po(r) , p = po(r) , v = Q(r)r , B = B (r) e z , (7) 

where (i = r, <p, z) are unit coordinate vectors, and we use the subscript to refer to the 
equilibrium. Equations (1), (3)-(6) are automatically satisfied by such an equilibrium, while 
equation (2) gives an additional requirement: 

1 dp t o 

9es = j- • 8 

Po dr 

Here, the effective gravitational acceleration g eS and the associated frequency fl e s are defined 
by 

g eS = fi e > = ^-n 2 r, (9) 
and the total pressure p t is given by 

Pt = P+^. (10) 

Equation (8) states that the net radial force, obtained by summing the gravitational, centrifu- 
gal and (total) pressure gradient forces, vanishes in equilibrium. Note that, to be consistent 
with the assumption of cylindrical symmetry, the background gravitational potential $ must 
be a function of only the cylindrical radius r. This is one of the many simplifications we 
have made in the model. 

As explained in §2.1, the initial equilibrium consists of two distinct zones separated at 
a radius r m (Fig. 1). At r m , there is a discontinuity in p and B, and sometimes also in fl 
and v (see §5 below). However, the pressure p is always continuous. 
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Linear Perturbations 



3.1. Wave Equation 



We now consider linear perturbations of the equilibrium configuration. We take the 
perturbations to be proportional to exp[i(—ut+m<f))], where u> is the mode frequency, which is 
in general complex, and m is the number of nodes in the direction; the azimuthal wavevector 
is equal to m/r. We assume that the perturbations are independent of z. Therefore, magnetic 
field lines which are initially parallel to z (by assumption) remain so even in the perturbed 
state, and the r- and ^-components of the magnetic field vanish at all times. We then have 
from equation (2), 



az 



. 



(11) 



Since we have further assumed that the gravitational potential $ is a function only of r, 
none of the force terms on the right-hand side of equation (2) has a vertical component. 
Therefore, the vertical velocity perturbation is also zero. The linear perturbations to the 
mass density, gas pressure, velocity, and magnetic field may be written as 



' Sf ' 
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Sp 


> = < 


' *v 


SB 





Pi(r) 
u(r) e r + v(r) e^ 
£i(r) e 2 



> x e 



—iuit+irruf) 



(12) 



Taking the curl of equation (2) and substituting equations (5) and (11) gives the following 
equation for the vorticity, 



d ,^ 1 
(V x v) = — Vp x 

at p 



v$ + 



dv 
~dt 



(13) 



Equations (5), (6), and (13) then form a complete set of equations for v and p. These 
equations are independent of the gas pressure and the magnetic field. The first order per- 
turbations of equations (5) and (6) give 

1 d . . im „ „. 

-— (ru) + — v = 0, 14 
r ar r 

dp0 n 

tapi — u —— = (J , (15) 
ar 

respectively, where 



a = lo — mO 



(16) 
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Similarly, the first order perturbation of equation (13) gives 

d ( PqK. 2 \ im 



Id. 

— (rpov) + 

rpo dr 



Pimg eS 

u = , 17 

Po ra 



apo dr \ 2Q J r 

where the epicyclic frequency k and the related vorticity frequency ( are defined by 

S4^ 2 w- < 18 > 

Finally, from equations (14), (15), and (17), a wave equation can be derived for the 
quantity W = ru: 



1 d ( dW\ 



m 



2 



rp dr \ ^° dr ) r 2 



tt 2 eS dlnpo 2( d\n(p () 



+ 1 



W = 0. (19) 



a 2 dlnr ma d\nr 

This is the fundamental equation for the model considered in this paper. The mode fre- 
quency uj is an eigenvalue of this second order differential equation and is determined by 
solving the equation and applying the boundary conditions. Any eigenvalue uj with a posi- 
tive imaginary part represents an unstable mode. We focus on such modes in the rest of the 
paper. It is possible to derive a necessary condition for the existence of instability, as shown 
in Appendix A, but we do not use that condition in the main paper. 

Once the eigenvalue u and the eigenfunction W = ru are obtained, we may solve for 
the other perturbed quantities by going back to the linear perturbation equations. Thus, 
equation (14) gives the perturbed azimuthal velocity v(r), and equation (15) gives the per- 
turbed density pi(r). To calculate the perturbed total pressure Pti(r), we use the azimuthal 
component of the first order perturbation of the momentum equation (2), 

2Qu — lav = . (20) 

r Po 

The perturbed magnetic field B x can be calculated from the first order perturbation of the 
induction equation (3), 

Bi = —u, 21 

a dr 

where equation (14) has been substituted. Finally, the perturbed gas pressure pi can be 
calculated from equations (20) and (21), by making use of p 1 = p tl — BqBi/A-k. 



3.2. Jump Conditions Across the Boundary 

Since the equilibrium model has a discontinuity at r = r m , where the accretion disk 
meets the magnetosphere, the solution must satisfy certain junction conditions at this bound- 
ary. Let us define r m+ = r m + + and r m _ = r m + 0~ as the radius r m on the two sides of 
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the boundary. Let us also define A m (/) to denote the jump in the value of any quantity / 
across the boundary, i.e. 

A m (/) = /(r m+ )-/(r m _) . (22) 

Consider first the Lagrangian displacement in the radial direction: £ — iu/a. Clearly, 
the displacement has to be the same on the two sides of the boundary. This gives the first 
jump condition: 

A m {^j = . (23) 

Since a is in general not continuous at r = r m (because Q need not be continuous), note 
that it is the quantity W/a that is continuous at the boundary, while W itself may or may 
not be continuous, depending on the behavior of Q. 

Consider next the total Lagrangian pressure perturbation at the displaced position of 
the fluid: 

Pti,L = Pa + £^ = Pa ~ £po9eS , (24) 

where we have used equation (8) in the last step. We may rewrite p t i in terms of u and v 
using equation (20) and then replace these in terms of W and dW/dr. Carrying out these 
steps, and requiring the Lagrangian pressure perturbation to be equal on the two sides of 
the boundary, we obtain the second jump condition: 



A m 



dW m ( mtt 2 cS 

crpo-; — - 2C PoW 

dr r V a 



. (25) 



The two jump conditions (23) and (25) form a complete set of relations that must be 
satisfied at r = r m . 



4. Models with Continuous Q at the Interface 

We first consider a model in which the angular velocity is continuous across r m , and 
take the profile of Q to be given by (see Fig. 2) 

Sl(r) = a m (—) 9 . (26) 



- 12 - 



The corresponding vorticity frequency is 



C(r) = (l - |) n(r) 



(27) 



For completeness, we consider values of the index q over the range from (constant angular 
velocity) to 2 (constant specific angular momentum). However, in real accretion disks, the 
relevant range is likely to be from 3/2 (Keplerian) to 2. We assume that the gravitational 
potential $ oc r _1 . 

A model with continuous Q is not likely to be relevant for accreting neutron stars, 
since we expect in general a discontinuity in Q at the interface between the disk and the 
magnetosphere. However, in the case of accretion onto a magnetosphere surrounding a black 
hole, since the field is not anchored to the black hole, we expect the magnetospheric fluid to 
corotate with the surrounding disk. Thus, the discussion in this section is most relevant for 
the black hole case. It is possible that the value of q may change across the boundary for a 
real black hole flow. We do not consider this possibility here, though the equations we write 
are general enough to handle it. 2 

As described in §2.1, we assume that there is a density jump at r rn . Thus, we model 
the density profile as 



where p + , p_, and 7 are constants. When 7 > 0, the mass density decays with increasing 
radius for r > r m . The jump in mass density at r m is described by the following dimensionless 
parameter, 



Note that /i can in general take any value between —1 and +1. However, we only consider 
values of /x > 0, i.e., p + > p_. 

Because Q is continuous at r m , so is a. Therefore, the two jump conditions written in 
§3.2 become simplified to 




(28) 



A m (W) = 



(30) 



2 Neglecting a possible change in q across the boundary is reasonable especially when the following fact 
is considered: when the disk has a sharp contrast in mass density across the boundary (i.e., p + ^> the 
solutions do not depend on the details of the disk inside r = r m . 
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dW 



dr 



m 2 n 2 



off 



2m( 



A m (p ) 



(31) 



The former equation states that W is continuous at r m . We write the value of W at this 
radius as W m . In the latter equation, p = p + for r = r m+ and p_ for r = r m _. 



4.1. Analytical Results 



We start first with the simplest case: (1) We assume that 7 = 0, i.e., the mass density 
is constant on the two sides of the boundary, with a jump at the junction described by the 
parameter p. (2) We assume that q is equal to either or 2. 

The first assumption implies that the term involving din p /dlnr in equation (19) van- 
ishes, while the second assumption means that the term involving <i ln(p C) l n r vanishes, 
either because ( is constant (for q = 0) or it is equal to (for q — 2). Therefore, away from 
the boundary r = r m , the wave equation (19) takes the form 



1 d 

r dr 



dr J 



m 



. 



(32) 



This differential equation has two simple solutions: W oc r +m , r~ m . We require the pertur- 
bations to decay away from the boundary, both as r — > oo and as r — > 0. Therefore, for 
m > 1, the solution must be of the form 



W = 



A _ r r, 



if r < r 



(33) 



A + r~ m , if r > r n 
where A + and A_ are constants. 

We now apply the two jump conditions (30) and (31). Since A m (W) = 0, we must have 

A + = A_r 2 ™ . (34) 

Thus, 

dW s 



(hr 



dr 



~(P+ + p-)mA_r. 



m—l 
m 



(35) 



Applying this to the second jump condition (31), we obtain a solution for the eigenvalue 
uj = a + mQ, 



mVL + p( ± yj ' —pmVL 2 cS + fi 2 ( 2 



9 = 0, 2, 



(36) 
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where \i is the density contrast parameter defined in equation (29). 



We are interested in unstable modes, i.e., modes with complex u>, where the growth 
rate of the instability is given by u\. From equation (36), we see that unstable modes exist 



whenever -/imfi^ + fi 2 ( 2 < 0, i.e. if 



mQ 2 „ 

<n<H! , A*i = — zf*- 



> . (37) 



c 2 

Thus, the instability exists only if ji is positive, i.e., if the density on the outside is greater 
than that on the inside. This is perfectly natural for the Rayleigh- Taylor instability. Sur- 
prisingly, if the density contrast is too large, i.e., if /i > fj,i, the instability shuts off. This 
is clearly the result of rotation, or more specifically vorticity. When q = 2 and the vorticity 
C = 0, then Hi — > oo, and the instability is present for any positive value of fi. However, for 
q = 0, we have ( = Q and /ii is finite. In this case, if f2^ ff small (i.e., the effective gravity is 
weak) and if we consider a low value of the azimuthal wavenumber m, the vorticity is able 
to eliminate the instability. 

The real part of oo gives the observed oscillation frequency of the mode. When q — 2, 
this is simply equal to mQ m . In this case, the mode is stationary in the frame of the fluid 
at the boundary, and what one observes is simply the Keplerian frequency Q m of the gas at 
the inner edge of the disk. However, when q = 0, the observed mode frequency is not equal 
to mQ m , but is equal to (m + fi)Q m , which differs from the orbital frequency. In this case, 
because of vorticity, we have a traveling mode in the fluid frame. 



4.2. Numerical Results 



For a more general disk, with 7 ^ and/or q 7^ 0,2, we have to solve the wave 
equation and the boundary conditions numerically. Although we derived in §3.1 a second- 
order differential equation as our basic wave equation (eq. [19]), for numerical purposes a 
pair of equivalent first-order differential equations is more convenient. Defining V = rv , we 
write 








1 



2r dC 
ma dr 



+ 



\ a 2 ma J 



K. ) L dpp 

ma J po dr 



1 dpi) 
Po dr 




(38) 



The first jump condition at r = r m is that W must be continuous (eq. [30]). From equa- 
tions (31) and (14), we have 

2C 
a 



A m (p o V0 



m n 2 cS 



W m A m p , 



(39) 
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where W m = W(r m ). Substituting equation (29) into this, we obtain the second jump 
condition 



[1 + H)V + - (1 - n)V- = 2i 



mQ 2 s 2C 



W m p . 



(40) 



We integrate the two first-order equations (38) from a small radius r = r\ <C r m outward, 
and simultaneously also from a large radius r = r 2 r m inward. Where the two integrations 
meet at r = r m , we adjust the relative normalization factor and the eigenvalue u such that 
the two jump conditions are satisfied. The starting values of W and V at r = r\ and r 2 
are chosen to correspond to the appropriate eigenvectors of the operator d/dr in equation 
(38). At r = ri, we require that the real part of the eigenvalue should be positive, so that 
d\W\ 2 /dr > and d\V\ 2 /dr > 0, while at r = r 2 , we require that the real part of the 
eigenvalue should be negative, so that d\W\ 2 /dr < and d\V\ 2 /dr < 0. The two eigenvalues 
of d/dr (see eq. [38]) are 



h,2 = 



1 dp 

Po dr 



1 dp 
po dr 
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ma dr 
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r dp 
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(41) 



The corresponding two eigenvectors (up to a normalization factor) are, respectively, 
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(42) 



It can be checked that, at both r = ri r m and r = r 2 ^> r m , k\ is positive and /c 2 is 
negative. Therefore, the eigenvector corresponding to k\ is a physical solution at r = r l5 and 
the eigenvector corresponding to fc 2 is a physical solution at r = r 2 . Thus, in our numerical 



calculations, we choose 

„-m+l 



and V = ir™ k\(r\)/m at r = ri, and 



and 



V = ir 2 m+ /c 2 (r 2 )/m at r = r 2 , as the starting values of W and V. 



One comment is in order regarding the integration. The point a = is a singularity 
of equation (38); it represents the corotation singularity. Therefore, while integrating the 
equation numerically, it is important to make sure that the integration contour goes above 
the singularity (Narayan, Goldreich, & Goodman 1987). When oj (thus a) is complex, this 
is easily arranged by just integrating along the real r-axis since the singularity is below the 
real axis for positive uj\ (which corresponds to an unstable mode, recall our convention that 
the perturbations are oc exp[— iut}). However, when uj is real, the corotation singularity is 
on the real r-axis. In this case, the integration contour must go above the singularity in the 
complex r-plane. 
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We present some numerical results in Figs. 3-6. In Figure 3, we show the frequency 
uj = u-R + iui, in units of Q m = Q(r m ), corresponding to a model with fx — 1, ^cfr,™ = 0-6^™, 
and 7 = (i.e., constant mass density away from r = r m ). The solid lines show uj\ and give 
the growth rate of unstable modes. The dashed lines represent cur, and correspond to the 
observed oscillation frequency of the mode. The results are perfectly reasonable — mode 
frequencies vary smoothly as a function of q between the two limits q = and q = 2 (where 
we have analytical results, §4.1). For m = 1 and small values of q, we see that the instability 
shuts off and we have instead two neutral modes with real uj. This is not unexpected. The 
analytical results in §4.1 indicate that for q = and this set of parameters, we have fx > fx\, 
which indicates that there should be no instability. 

For the parameters used in Fig. 3 and for q in the range 1.5-2 of interest for accretion 
disks, we find that all azimuthal wave numbers m are unstable. However, as equation (37) 
shows, whether or not a given mode is stable or unstable depends on whether fx is greater 
than or less than fX\ and this depends on the magnitude of VL 2 S /C, 2 . (Although the result 
given in eq. [37] is true only for q — 0, 2, the qualitative result is valid also for other values 
of q.) In Figure 4 we show the numerically determined critical values of Q^ s m that separate 
stable and unstable modes in the q — f^ ff m plane (we restrict our attention to q in the range 
1.5-2 which is of most interest for QPOs). The results correspond to fx = 1, 7 = 0. We 
see that, for q < 2, there always exists a range of Q 2 Sm over which modes with a given m 
are stable. For example, when q = 1.5 and 0.024 < Ql Sm /Q 2 n < 0.04, modes with m — 1 
are stable and modes with m > 2 are unstable, while for 0.018 < Vt 2 eS m /Vt 2 m < 0.024, modes 
with m < 2 are stable and modes with m > 3 are unstable. 

Figure 5 shows results for some other choices of parameters: fx = 0.8, fij^j m = an d 
7 = in panel (a); fi = 1, Q 2 eS:Jn = n 2 m , and 7 = in panel (b); and, fi = 1, n 2 cS m = tt^, and 
7 = 1 in panel (c). Comparing panel (b) with Figure 3, we see that the zone of stable modes 
at small q and m = 1 seen in the latter, disappears when fi^ ff increases. This is because \i\ 
increases and becomes larger than fi. All other features are very similar. Comparing panels 
(a) and (b) of Figure 5, we see the effect of varying fx. For this particular set of parameters, 
for m — 1 we see that uj\ decreases with increasing fx when q is close to 0, but increases with 
increasing fx when q is close to 2. For m > 1, however, ui\ always increases with increasing 
fx. The real part of the frequency, cur, increases with increasing fx for < q < 2, but remains 
independent of fx for q = 2 (see eq. [36]). Finally, by comparing panels (b) and (c), we see the 
effect of having a non-constant density (7 7^ 0). Generally, as 7 increases, uj\ also increases, 
but wr is hardly affected. 

In Figure 6 we show results corresponding to a Keplerian-type disk (q = 3/2), as a 
function of fx, for ^l 2 cSm = 0.2f2^ and 7 = 1 (i.e., the mass density decays with increasing 
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radius according to p oc r _1 ). We see that the imaginary part of the frequency uj\ strongly 
depends on /i, with u\ — > as fi — > 0. This is because the instability is driven by the density 
contrast at the boundary. As the density contrast goes to zero, the instability must clearly 
vanish. The real part of the frequency, wr, however, depends only weakly on \i. 

Summarizing, we find that the Rayleigh- Taylor instability operates over a wide range 
of parameters so long as the density on the outside is greater than the density on the inside. 
The growth rate of the instability is typically large, of order ~ f2 e g, which means that the 
growth occurs fairly rapidly. For 3/2 < q < 2, the range of interest for accretion disks, the 
real part of the frequency is nearly equal to mQ m in almost all cases. This is because the 
mode nearly corotates with the gas at r m . For a large density contrast // — > 1 and a weak 
effective gravity fll s <C Q^, low-m modes become stable and the instability is limited to 
higher-m modes. 



We now consider a model that is more appropriate for an accretion flow around a 
magnetized neutron star (Ghosh & Lamb 1978, 1979). Since the magnetic field of a neutron 
star is frozen to the star, the low density plasma in the magnetosphere (r < r m ) must corotate 
rigidly with the star. However, outside the magnetospheric radius, we expect the accreting 
gas to rotate differentially as in a standard accretion disk. Also, the angular velocity of 
the disk at r m will, in general, not match the angular velocity of the star. Based on these 
considerations, we assume the following model for the angular velocity (Fig. 7): 



where Q + and Q_ are constants. In general Q_ ^ Q + . For the density profile, we assume 
the same model as in equation (28). 

The jump conditions are given by equations (23) and (25). 



As in §4.1, we consider the special case when 7 = and q = or 2. We may then 
repeat the same steps described in §4.1, except that the jump conditions are now different. 
We obtain the following results for the eigenvalue, 



5. 



Models with Discontinuous f2 at the Interface 




(43) 



5.1. Analytical Results 



U = 77 



(1 + fj) (mn + + C+) + (1 - /i) (mtt_ - (-) ± 



(44) 
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where the "+" and the "—"in the subscripts denote evaluations at r = r m+ and at r = r, 
respectively, and 



When VL + = fi_ (which means ( + = ^cfr,+ = ^eff,-), equation (44) simplifies to equa- 
tion (36). 

When 6 is positive, there is an instability. A study of the expression for 6 indicates that 
there are two distinct terms, written in separate lines in equation (45). Therefore, there are 
two distinct mechanisms of instability. The first term in (45) is the one we have already seen 
in §4, due to the Rayleigh- Taylor instability. If Q and its derivative are continuous at r m , this 
is the only term present, and it gives an instability only if p + > p_, i.e., /i > 1. However, for 
the problem at hand, Q is discontinuous at r m , and so the second term can also be important. 
Ignoring the vorticity terms (± for simplicity, this term is directly proportional to (Q + — Q ) 2 , 
i.e., the square of the angular velocity jump across the boundary. This term is nothing but 
the classical Kelvin-Helmholtz instability associated with a velocity discontinuity. 

The observed mode frequency is given by the real part of uo in equation (44). If p = 1, 
i.e., we have the maximum density contrast across the boundary, then the frequency is just 
equal to the appropriate frequency mVt + + ( + of the outer accretion disk at r m+ . That is, 
the rotation rate of the neutron star is irrelevant. This is reasonable since there is no 
mass associated with the flow inside r m , and so the frequency f2_ of this region of the flow 
should have no influence on the mode. However, when /i ^ 1, the mode frequency is a linear 
combination of mVt + + ( + and mfl- + with weights given by the two densities. 

All of these results are very reminiscent of the classical results for the Rayleigh- Taylor 
and Kelvin-Helmholtz instability (Drazin & Reid 1981). We discuss this connection in Ap- 
pendix B. 



For the general case with 7 ^ and/or q ^ 0,2, we numerically solve equations (38), 
(23), and (46), following the numerical approach described in §4.1. We write the jump 
condition (25) as 



= 2[(l + /U )(mfi c 2 ffi+ -a)-(l-/i)MeV + C-)] 



(45) 



5.2. 



Numerical Results 




(46) 
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In Figure 8 we show some numerical results for a disk with 7 = 0, [i — 0.4 and the 
frequency of effective gravity fll s + = In panel (a), the disk angular velocity satisfies 
f2_ = 2Q + , i.e., the neutron star rotates twice as fast as the disk at r = r m . In panel (b), 
the disk angular velocity satisfies f2_ = Q+, i.e., the neutron star rotates at the same rate 
as the disk and the angular velocity is continuous (but not smooth) at r = r m . In panel 
(c), the disk angular velocity satisfies = 0.4f2 + , i.e., the neutron star rotates more slowly 
than the disk. 

Comparing the three panels, we see that the real part of the frequency (wr) monotoni- 
cally decreases as decreases. This is as expected. Based on the discussion in §5.1, the 

mode frequency is a linear combination of the outer and inner frequencies at the boundary. 
As the latter decreases, the mode frequency must also decrease. The growth rate of the mode 
is more interesting; it is minimum when = 1, and increases both for fi_/fi + > 1 and 

fi_/fi+ < 1. This behavior is because of the simultaneous presence of the Rayleigh- Taylor 
and Kelvin-Helmholtz instabilities. When = 1, there is no velocity discontinuity at 

the boundary and we have only the Rayleigh- Taylor instability, with a certain growth rate. 
However, both for Q-/Q + < 1 and > 1, there is an additional contribution to the 

growth rate from the Kelvin-Helmholtz instability, and so the net growth rate is enhanced. 

Note that, for the case shown in panel (a), we have f^ ff _ = — 0.5f2?_, indicating that 
the disk is super-Keplerian right inside r = r m . But there are still unstable modes in this 
case, since the disk is sub-Keplerian for r > r m . For the cases shown in panels (b) and (c), 
the disk is sub-Keplerian for both r < r m and r > r m . 

6. Summary and Discussion 

The main message of this paper is that there are instabilities associated with the mag- 
netospheric radius where an accretion disk meets the magnetosphere of the central mass and 
that these instabilities may be relevant for understanding QPOs in binary systems. There 
are two natural instabilities at the disk-magnetosphere interface: (1) Rayleigh- Taylor (or 
interchange) instability associated with a density jump, and (2) Kelvin-Helmholtz instabil- 
ity associated with an angular velocity jump. These instabilities, which are well-known for 
classical uniform Cartesian flows (Appendix B), survive with relatively little modification for 
a rotating, shearing flow with a density gradient (§§4,5). The unstable modes are expected 
to grow to become nonaxisymmetric perturbations with large amplitude (e.g., the streams 
in the simulations of Igumenshchev et al. 2003) and to give strong quasi-periodic variations 
in the observed intensity. 
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If, as we propose, these instabilities are responsible for the observed QPOs in X-ray 
binaries, the same model might work both for neutron stars and black holes. The existence 
of a disk-magnetosphere interface around accreting magnetized neutron stars is well-known. 
The popular beat frequency model for QPOs invokes blobs in the accretion disk orbiting at 
the Keplerian frequency at the magnetospheric radius (Alpar & Shaham 1985; Lamb et al. 
1985; Strohmayer et al. 1996) or at the sonic radius (Miller, Lamb, & Psaltis 1998; Lamb 
& Miller 2001). While the model does not explain the origin of the blobs, it is reasonable 
to assume that the blobs are in some cases at least created by the interchange and Kelvin- 
Helmholtz instabilities studied in this paper. A fact not widely appreciated is that the 
magnetospheric model might also apply to accreting black holes and unmagnetized neutron 
stars. As Bisnovatyi-Kogan & Ruzmaikin (1976) showed, it is possible for the region close to 
a black hole to become magnetically dominant (Livio et al. 2003; Narayan et al. 2003). An 
accretion disk would then be disrupted at a magnetospheric radius, just as in the neutron 
star case, and the interface would be unstable and produce QPOs. The main difference 
between the two cases is that the magnetic field is not anchored to the black hole, and so 
the magnetosphere does not rotate rigidly with the star as in the neutron star case. This 
difference between the black hole and neutron star problems leads to some differences in the 
results for the two cases, as discussed in §§4,5. 

Our analysis shows that nearly all azimuthal wavenumbers m are unstable under rea- 
sonable conditions. Although modes with higher values of m grow more rapidly, we expect 
that these modes will saturate at relatively small amplitudes. The low-m modes, on the 
other hand, are likely to grow to large amplitude and are therefore of most interest for un- 
derstanding QPOs. For q in the range 3/2 to 2 (Keplerian to constant angular momentum), 
and reasonable assumptions about the density profile, we have shown that the observed 
mode frequency tends to be of order mO m . Thus, in the simplest version of the model, we 
expect the observed QPO frequencies to be in the ratio 1:2:3, etc. However, as we showed 
in §4 (see Fig. 4), it may often be the case that the m — 1 mode is stable and that only 
modes with m > 2 are unstable. This should happen whenever the effective gravity in the 
radial direction is weak, i.e., when the gas pressure in the disk is low. In this case, the QPO 
frequencies should be roughly in the ratio 2:3:4, etc. It is interesting that a frequency ratio 
2:3 is frequently seen in both black hole and neutron star systems (Abramowicz & Kluzniak 
2001; Abramowicz et al. 2003; Remillard et al. 2002a). 

We have emphasized that, in our opinion, QPOs should be produced by the growing 
modes in the disk. Our belief is based on the fact that the amplitude of intensity fluctuations 
observed in QPOs is often quite large. Observations also show that QPOs have a finite 
frequency width, which indicates that the unstable modes in the disk have a finite life time 
(hence the term quasi-periodic oscillations). Apart from the rotation period, there are two 
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natural time scales in the disk: the time scale associated with the effective gravity, which is 
~ l/ficff) an d the viscous time scale, which is ~ r/v r , where v r is the mean radial velocity 
of the disk fluid. For a standard disk, Q c g ~ c s /r, where c s is the sound speed. Since 
c s ^> v r , the time scale l/fi c fr is generally much shorter than r/v r . Therefore, the likely 
lifetime of blobs created by the instability is 1/Q c g. Once a mass blob is formed, it will drift 
toward the central object under the action of gravity, causing a displacement in the center 
frequency and a width to the QPO feature in the power spectrum. The width is likely to 
be Af ~ il cS /2n ~ (/i/r)(Q m /27r), where h is the vertical thickness of the disk, and the 
displacement speed is approximately df/dt ~ (1^/271-)/. 

Another issue concerns how the presence of a nonaxisymmetric mode translates to a 
time modulation of the observed flux. Two possibilities are likely. One is that the system 
is viewed in a nearly edge-on configuration so that the accreting star eclipses the far side of 
the disk. Then, as bright and faint segments of the disk are successively eclipsed the signal 
at the observer will be modulated. The other possibility, which also requires fairly high 
inclination, is that the motion of the gas is relativistic and the observed signal is dominated 
by the blue-shifted segment of the disk. Once again, as the nonaxisymmetric pattern rotates, 
the signal will oscillate. Both mechanisms require that the bright and faint patches on the 
disk should have large areas, since otherwise the fractional modulation of the observed X-ray 
flux will be small. 

The azimuthal extent of a bright patch in a nonaxisymmetric mode with wavenumber m 
is ~ 7r/m. The radial extent also has an m-dependence, since away from the magnetospheric 
radius the perturbation solutions decay with radius as ~ r ±m (§4.1). For both reasons, low-m 
modes have patches with the largest area and hence are most promising. The area occupied 
by a bright patch, defined to correspond to an annular region in the disk bounded by a radius 
where the perturbation amplitude is half the peak, is estimated to be S ~ (7rr^/2m)(4 1//m — 
4 -1 / m ). For m= 1,2,3 we have S ~ 1.8757rr^, 0.3757rr^, 0.167rr^, respectively. For large m, 
S approaches zero according to S ~ 27rr^(ln2/m 2 ). The scaling clearly shows that low-m 
modes dominate by a large factor. In addition, as we argued earlier, low-m modes are likely 
to saturate with substantially larger amplitudes than high-m modes. This is yet another 
reason why only the lowest order few modes are expected to cause a discernible signal in the 
observations. 

Apart from demonstrating that unstable modes exist and that mode frequencies in the 
ratio 2:3 are possible, the model does not really explain any of the many puzzling features seen 
in the observations (see the summary in §1). For instance, the model does not explain why 
the frequency difference between twin kHz QPOs in neutron star systems is often roughly 
of order the neutron star spin frequency (van der Klis 2000) or sometimes half the spin 
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frequency (Wijnands et al. 2003). Even though the version of the model described in §5 
appears to have the necessary ingredients for the beat frequency model to operate, namely gas 
orbiting at Keplerian frequency around a magnetosphere that rotates at the stellar frequency, 
nevertheless our analysis does not reveal any beat phenomenon. Clearly, additional physics 
is needed beyond what we have considered here. 

A Rayleigh- Taylor- like process has been studied extensively by Titarchuk and collabo- 
rators in a sub-Keplerian transition region of a disk around a black hole or a weakly mag- 
netized neutron star (Titarchuk, Lapidus, & Muslimov 1998; Osherovich & Titarchuk 1999; 
Titarchuk, Osherovich, & Kuznetsov 1999; Titarchuk 2002, 2003). These authors focus only 
on stable modes and suggest that their model can explain the observed correlation between 
the twin kilohertz frequencies and the horizontal branch QPO frequency (Osherovich & 
Titarchuk 1999; Titarchuk 2003). In their model, the dynamical effect of the magnetic field 
is always assumed to be unimportant, so the fluid is described purely within hydrodynamics. 
The model assumes the existence of a thin sub-Keplerian transition region in the vicinity 
of the compact central object where the accreting matter adjusts itself either to the surface 
of a rotating neutron star or to the innermost boundary of the accretion disk (Titarchuk, 
Lapidus, & Muslimov 1998). But the origin of the transition layer is not explained, especially 
considering that the magnetic field is assumed to be weak. 

The model described in the present paper (see §2.1) differs from Titarchuk's model in 
two respects. First, we assume that the magnetic field is dynamically important inside the 
inner edge of the disk. Second, we focus on genuinely unstable modes rather than on stable 
modes, since we assume that only unstable modes can grow to a large enough amplitude to 
produce the observed intensity fluctuations. Because of the assumption of a strong magnetic 
field, a narrow transition region develops naturally at the disk-magnetosphere interface of 
the model, and unstable oscillations are triggered at this interface. 

Kaisig et al. (1992) used two-dimensional shearing box simulations to study the magnetic 
interchange instability in an accretion disk with a vertical magnetic field. When the magnetic 
field is strong and its strength decreases with increasing radius, they find that an instability 
develops spontaneously. The model we consider is an extreme version of the Kaisig et al. 
model in which the magnetic field decreases discontinuously at the magnetospheric radius. 
Our analytical results are consistent with Kaisig et al. (1992) 's numerical result that, in 
the linear regime, the growth rate of perturbations depends on the azimuthal wave number 
k = m/r of the initial perturbations as uj\ ~ k 1 ^ 2 . In both studies, the dependence of 
quantities on z is neglected. When this dependence is included, Lovelace, Romanova, & 
Newman (1994) showed that the twisting of field lines acts to drive winds or jets from the 
disk surfaces, which increases the disk accretion speed and so amplifies the magnetic field and 
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leads to runaway or implosive accretion and explosive wind or jet formation. Similar results 
have also been obtained by Lubow, Papaloizou & Pringle (1994a,b). [However, Lubow et al. 
1994b also claimed the existence of stable solutions with no accretion and no wind.] 

Lubow & Spruit (1995) and Spruit et al. (1995) found that disk shear tends to stabilize 
the disk-magnetosphere interface. This seems to be in conflict with our result that the growth 
rate of the Rayleigh- Taylor instability increases with increasing q. However, we stress that 
in their model the magnetic field has both vertical and radial components. Indeed, they 
modeled the disk as a sheet of zero thickness, so that the radial component of the magnetic 
field has a jump as one goes from below to above the disk. The pressure of gas and magnetic 
fields is negligible in the disk, but the magnetic curvature force appears in the equation 
of motion. This is very different from our model, where the disk has a nonzero (indeed 
infinite) thickness, the magnetic field has only a vertical component, the curvature force of 
the magnetic field is zero, and the magnetic pressure and gas pressure both play a dynamical 
role. Our model is perhaps more suitable for describing the central layer of the disk, while 
their model may be more applicable to the surface layers. 

Chandran (2001) has considered the Rayleigh- Taylor instability of a strong vertical mag- 
netic field confined by a disk threaded with a horizontal magnetic field. The aim of his study 
was to understand the equilibrium of magnetic flux tubes observed in the central regions of 
the Galaxy. Although the model has some points of similarity with the present study, there 
are also large differences. The disk in Chandran's model is assumed to have a uniform rota- 
tion and the gravitational potential is assumed to correspond to a constant background mass 
density (so that the gravitational acceleration increases linearly with radius). In contrast, 
we assume that the disk is differentially rotating and that the gravitational potential corre- 
sponds to that of a compact mass at the center. However, Chandran considers a compressible 
gas whereas we simplify our problem by taking the gas to be incompressible. 

Chandran (2001) has derived a formula (his eqs. [91] and [92]) for the oscillation fre- 
quency when there is no magnetic field and the density contrast parameter \i = ±1. His 
results are consistent with our analytical results for a disk with constant angular velocity 
(see our eq. [36]). In particular, the results confirm that the disk vorticity has the effect of 
stabilizing the modes. 
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A. A Necessary Condition for Disk Instability 



A necessary condition for instability can be derived from the wave equation (19), fol- 
lowing the approach pioneered by Howard (1961) (see also Li et al. 2003, Appendix B). 



Defining W = ipyfa, equation (19) can be recast as 
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whole flow yields 



rr2 
J ri 



dip 




dr 

2C ^ln(poC) 
m dlnr 



dtt\ 2 ' 
dr~) 



1 

- + a 
a 



2mp dr 



TPo ~~d~r) ^ ^ ' TP ° dr 



*d<p 
rpocrip — 
dr 



r-2 



n 



(A2) 



Taking the boundary condition to be W = at r\ and r2, we see that the right-hand side 
of equation (A2) vanishes. Then, the imaginary part of the left-hand side of equation (A2) 
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where cuj = $s(u) = ^(c). The integral in equation (A3) is positive, and therefore uj\ = 0, 
unless the condition, 

1 



4 ' 



is satisfied somewhere in the region of integration, where 

(I 111 /)n / (l\ ' \ ' " ' 
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din p 
dlnr 
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cS 



dlnfl 
dlnr 



din po 
dlnr 



(A4) 



(A5) 



Note that f2g ff can be either positive or negative. 



Therefore, a necessary condition for the existence of instability (i.e., u\ ^ 0) is that the 
inequality (A4) must hold somewhere in the disk. Obviously, the converse of the inequality, 
i.e. R e s > 1/4, is a sufficient condition for stability. 



3 If the boundary condition is taken to be dW/dr = 0, then it can be shown that the right-hand side of 
eq. (A2) is real, and the result is still the same. 



-25 - 



B. Classical Ray leigh- Taylor and Kelvin-Helmholtz Instability 



Let us start from the hydro dynamic equations (1) and (2) (with B = 0), and use Carte- 
sian coordinates (x,y,z). If the fluid is incompressible, equation (1) becomes equations (5) 
and (6). Let us assume that the system is two-dimensional and that all quantities are inde- 
pendent of z. We take gravity to be in the x direction: g = — V$(rr) = —g(x)e x . For the 
unperturbed equilibrium state we assume that p = po(x), p = po(x), and v = V(x)e y . Then, 
the equilibrium in the x-direction is given by 



9 



1 dp 
p dx 



(Bl) 



Consider now linear perturbations 




pi 0*0 

Pi(x) > x f 

u(x) e x + v(x) e y 



■iwt+iky-y 



(B2) 



where u> and k y are constants. The first order perturbations of equations (5) and (6) then 
give 



respectively, where 



du ., 

dx y 
dpo 

la — u— — 

dx 



a = uj — k y V 



0, 
0, 



(B3) 
(B4) 

(B5) 



The first order perturbation of equation (13), which is derived from equation (5) and the 
curl of equation (2), leads to 



1 d 

Po dx 



(Pov) + 



id/ dV\ 
apo ax \ ax 



pi k y 
u = —- -g 

Po cr 



(B6) 



^From equations (B3)-(B6) we can derive a wave equation for u, 

g dlnp 2( d\n(p () 



1 d ( du\ 2 
J«dx\ P »Tx)- k * 

where the vorticity ( is given by 



a 2 x dlnx 



kyXa dlnx 
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As in the main paper, we can derive the following two junction conditions at a surface 
of discontinuity at x = 0: 

A m (£) = , (B9) 



du , / k v g 
ap0 dx-- k «{ — 



2( p u 



(BIO) 



where A m (/) = f(x = 0+) - f(x = 0"). 



Following Drazin & Reid (1981, Chapter 4, except that we switch x and z to be consistent 
with the convention in this paper), we assume for the unperturbed equilibrium state 



P- , 
P+ , 



if x < 
if x > 



v = 



V_e y , if x < 
y + e y , if x > 



(BH) 



where p_, p + , and are constants. Then, on each side of x = 0, we have £ = and 
dp/dx = 0, and the solutions of the wave equation (B7) are u oc exp(±k y x). Requiring that 
u(x — > ±oo) = 0, then from the junction conditions (B9) and (B10) we can solve for the 
eigenvalue for the frequency (assuming that k y > 0) to obtain 



uj = k. 



p + V+ + P-V^ 
P+ + P- 



± 



P+-P- 
P+ + P- 



k y g 



P+P- 
P+ + P- 



kl{V + -V-f 



(B12) 



Equation (B12) is exactly the same as equation (4.20) in Drazin & Reid (1981), provided 
we substitute their s with our — icu, and their subscripts "1" and "2" with our subscripts 
"— " and "+", respectively. The first term under the square root in equation (B12) repre- 
sents the Rayleigh- Taylor instability, while the second term represents the Kelvin-Helmholtz 
instability. 

The analysis presented in the main paper goes beyond this classical analysis in a few 
respects. We include rotation (and correspondingly vorticity), our flow has shear (actually 
differential rotation), and we consider a density gradient. These additional features modify 
some of the results, but the basic physics remains the same, namely the presence of the 
Rayleigh- Taylor and Kelvin-Helmholtz instabilities. 
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Fig. I. — Schematic representation of the (r, z)-section of the cylindrical disk model studied 
in the paper. The wavy line represents the central compact star. The two thick vertical 
lines represent the radius r = r m which separates the magnetosphere on the inside from the 
accretion disk on the outside. The thin vertical lines with arrows show magnetic field lines. 
In the region r > r m , the mass density (p + ) is high but the magnetic field is weak. In the 
region r < r m , the mass density is low but the magnetic field is strong. The jump in 
the mass density at r = r m is described by the parameter ji defined in eq. (29). 
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Fig. 2. — Disk angular velocity as a function of radius for the case of a central black hole 
(§4.1). The vertical dashed line marks the magnetospheric radius r m . The angular velocity 
is taken to be continuous across this boundary. 
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Fig. 3. — Numerical solutions for the mode frequency as a function of q (defined by f2 oc r~ q ). 
The mass density has a jump at r = r m corresponding to /i — 1, and the density is assumed to 
be constant on the two sides of the boundary. The effective gravity is given by Qg ff m = 0.6f2^. 
Solid lines correspond to the imaginary part of the frequency and measure the growth rates 
of the modes. Dashed lines correspond to the real part of the frequency and represent the 
observed oscillation frequencies. The different lines correspond to m — 1, 2, 3, 4, and 5 
(bottom to top). 
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1.5 1.6 1.7 1.8 1.9 2 

q 

Fig. 4. — Critical effective gravity frequency ^e ffm , which separates stable and unstable 
modes, for 3/2 < q < 2. It is assumed that /i = 1 and 7 = 0. The different lines correspond 
to m — 1, 2, 3, 4, and 5 (top to bottom). Each line separates the q-ill Sm plane into two 
regions: above the line modes of the given m are unstable, and below the line the modes are 
stable. 
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(b) 




(c) 



Fig. 5. — Similar to Fig. 3, but (a) fi = 0.8, ^cffm = an d the mass density is constant 
on both sides of r = r s ; (b) /i — 1, fi;? ffm = Q 2 rn1 and the mass density is constant on both 
sides of r = r m ; (c) /z = 1, ^effm. = an d the mass density is constant for r < r m , but 
decays with increasing radius according to p oc r _1 for r > r m . 
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Fig. 6. — Mode frequencies as a function of \x (defined by the jump in mass density at r = r m ) 
for a Keplerian-type disk (q = 3/2). The mass density is constant for r < r m , but decays 
with increasing radius according to p oc r _1 for r > r m . At r = r rn the effective gravity 
frequency is Q 2 cSm = 0.2Q^. Solid lines correspond to the imaginary part of the frequency 
and dashed lines to the real part of the frequency. The different lines correspond to m — 1, 
2, 3, 4, and 5 (bottom to top). 
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Fig. 7. — Disk angular velocity as a function of radius for the case of a magnetized neutron 
star (§4.2). The vertical dashed line marks the magnetospheric radius r m . The angular 
velocity has a jump at r rn : A m (f2) = Q + — f2_ ^ 0. The jump may be either positive or 
negative. 
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Fig. 8. — Mode frequencies for a disk with an angular frequency jump at T m ~. Q — — 
constant for r < r m , and Q = Q + (r/r m )~ q for r > r rn . The mass density is constant away 
from r = r m , but has a jump at r m with /i = 0.4. At r = r m+ , it is assumed that 



eff,+ 



+ ■ 



(a) The angular velocity is assumed to satisfy f2_ = 2f2 + . (b) = (c) = 0.4f2^ 
Line types have the same meanings as in Fig. 3. 



